clear,clc,clf
b=[1,-3,13,27,18];
a=[15,11,2,-4,-2,1];
N=30;
delta=impz(b,a,N);
x=[ones(1,5),zeros(1,N-5)];
h=filter(b,a,delta);
y=filter(b,a,x);
subplot(121);stem(h);title('h(n)')
subplot(122);stem(y);title('y(n)')
%% 7.3.1 巴特沃斯滤波器
% 1. batter函数
% batter函数用于计算N阶巴特沃斯归一化（3dB截止频率omega c=1）模拟低通原型滤波
% 器系统函数的极点和增益因子
% [z,p,k]=buttap(N)
% 返回极点和增益
% z=一个空矩阵
% p=设计的巴特沃斯滤波器的零点
% k=设计的巴特沃斯滤波器的增益
% N=设计的巴特沃斯滤波器阶数

% 2.buttord函数
% [m,Wn]=buttord(Wp,Ws,Rp,Rs)
% Wn：3dB截止频率；Wp：通带截止频率；Ws:阻带起始频率；
% Rp：通带内波动；Rs：阻带内最小衰减

% 3. butter函数
% [b,a]=butter(n,Wn)
% 返回巴特沃斯低通滤波器的系数
% b:H(z)分子多项式系数
% a:H(z)分母多项式系数
% [b,a]=butter(n,Wn,'high')
% 返回巴特沃斯高通滤波器的系数
% [b,a]=butter(n,[W1,W2])

%% 7-12
clear,clc,clf
[z,p,k]=buttap(20);
[num,den]=zp2tf(z,p,k);
freqz(num,den);

%% 7-13
w1=85/500;
w2=125/500;
[B,A]=butter(1,[w1,w2],'stop');
[h,w]=freqz(B,A);
f=w/pi*500;
plot(f,20*log(abs(h)));
grid
axis([50,150,-30,10]);
xlabel('频率/Hz');ylabel('幅度/dB')

%% 7-14
% 设计模拟巴特沃斯低通滤波器，并绘制幅频特性响应曲线
clear,clc,clf
n=0:0.01:2;
for i=1:4
    switch i
        case 1
            N=1;
        case 2
            N=3;
        case 3
            N=8;
        case 4
            N=12;
    end
    [z,p,k]=buttap(N);
    [b,a]=zp2tf(z,p,k);
    [h,w]=freqs(b,a,n);
    magh=abs(h);
    subplot(2,2,i);plot(w,magh)
    axis([0 2 0 1]);
    xlabel('w/wc');ylabel('|H(jw)^2');
    title([' filter N = ',num2str(N)]);
    grid on
end
%% 7-15
clear,clc,clf
Fs=40000;fp=5000;fs=9000;
rp=1;rs=30;
wp=2*fp/Fs;ws=2*fs/Fs;
[N,wc]=buttord(wp,ws,rp,rs);     % 计算数字滤波器的阶数和通带频率
[b,a]=butter(N,wc);
w=0:0.01*pi:pi;
[h,w]=freqz(b,a,w);
plot(w/pi,20*log10(abs(h)),'k');
axis([0,1,-100,10]);grid
xlabel('\omega/\pi');ylabel('幅度/dB');
%% 设计一个低通滤波器，采样速率为10kHz，fp=2kHz,fs=3kHz,Rp=4dB,Rs=30dB
% 7-16
clear,clc,clf
fn=10000;fp=2000;fs=3000;
Wp=fp/(fn/2);
Ws=fs/(fn/2);
Rp=4;
Rs=30;
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
[H,F]=freqz(b,a,1000,8000);% freqz(b,a,计算点数，采样速率)
subplot(121);plot(F,20*log10(abs(H)));grid on
xlabel('频率/Hz');ylabel('幅值/dB');
axis([0 4000 -30 3]);
pha=angle(H)*180/pi;
subplot(122);plot(F,pha);grid on 
xlabel('频率/Hz');ylabel('相位');

%% 设计一个高通滤波器，采样速率为10kHz，fp=0.9kHz,fs=0.6kHz,Rp=3dB,Rs=20dB
% 7-17
clear,clc,clf
fn=10000;fp=900;fs=600;
Wp=fp/(fn/2);
Ws=fs/(fn/2);
Rp=3;
Rs=20;
[n,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn,'high');
[H,F]=freqz(b,a,900,10000);% freqz(b,a,计算点数，采样速率)
subplot(121);plot(F,20*log10(abs(H)));grid on
xlabel('频率/Hz');ylabel('幅值/dB');
axis([0 4000 -30 3]);
pha=angle(H)*180/pi;
subplot(122);plot(F,pha);grid on 
xlabel('频率/Hz');ylabel('相位');

%% 7.3.2 切比雪夫I型滤波器
% 1.cheb1ap
% cheb1ord
% cheby1
% 7-20 设计切比雪夫I型低通滤波器示例
clear,clc,clf
Wp=3*pi*4*12^3;
Ws=3*pi*12*10^3;
rp=1;
rs=30; % 设计滤波器的参数
wp=1;ws=Ws/Wp; %对参数归一化
[N,wc]=cheb1ord(wp,ws,rp,rs,'s'); % 计算滤波器阶数和阻带起始频率
[z,p,k]=cheb1ap(N,rs); % 计算零点、极点、增益
[B,A]=zp2tf(z,p,k); % 计算系统函数的多项式
w=0:0.02*pi:pi;
[h,w]=freqs(B,A,w);
plot(w*wc/wp,20*log10(abs(h)),'k');grid 
xlabel('\lambda');ylabel('A(\lambda)/dB')

%% 7-21
% 设计一个数字高通滤波器，它的的通带范围为400-500Hz，通带内容许有0.5dB的波动，
% 在小于317Hz的频带内阻带内衰减至少为19dB，采样频率为100Hz。
clear,clc,clf
wc=2*1000*tan(2*pi*400/(2*1000));
wt=2*1000*tan(2*pi*317/(2*1000));
[N,wn]=cheb1ord(wc,wt,0.5,19,'s');
[B,A]=cheby1(N,0.5,wn,'high','s');
[num,den]=bilinear(B,A,1000);%数字滤波器设计
[h,w]=freqz(num,den);
f=w/pi*500;
plot(f,20*log10(abs(h)));grid
axis([0,500,-80,10])
xlabel('频率/Hz');
ylabel('幅度/dB');

%% 7-22
clear,clc,clf
Wp=[60 200]/500;
Ws=[50 250]/500;
Rp=3;
Rs=40;
[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
[b,a]=butter(n,Wn);
freqz(b,a,128,1000);
title('n=7 巴特沃斯滤波器')

%% 7-23
clear,clc,clf
n=0:0.02:4;
for i=1:4 % 取4种滤波器
    switch i
        case 1; N=1;
        case 2; N=3;
        case 3; N=5;
        case 4; N=7;
    end
    Rp=1;%通滤波纹为1dB
    [z,p,k]=cheb1ap(N,Rp);
    [b,a]=zp2tf(z,p,k); %将零点极点增益形式转换为传递函数形式
    [H,w]=freqs(b,a,n); %按n指定的频率点给出频率响应
    magH2=(abs(H)).^2;    %给出传递函数幅度平方
    subplot(2,2,i);
    plot(w,magH2);grid on
    title(['N=' num2str(N)]);
    xlabel('w/wc');ylabel('H(jw)|^2')
end

%% 7-24
% 设计一个6阶切比雪夫II型模拟低通滤波器，阻带波动为70dB，并显示其幅值和相位响应
clear,clc,clf
[z,p,k]=cheb2ap(6,70);  %低通滤波器原型
[num,den]=zp2tf(z,p,k);
freqs(num,den)
%% 7-25
clear,clc,clf
Wp=40/500;
Ws=150/500;
Rp=3;
Rs=60;
[n,Ws]=cheb2ord(Wp,Ws,Rp,Rs)
[b,a]=cheby2(n,Rs,Ws)
freqz(b,a,512,1000)
title('切比雪夫II型低通滤波器')

%% 7-26
clear,clc,clf
Wp=3*pi*4*12^3;
Ws=3*pi*12*10^3;
rp=1;
rs=30;
wp=1;ws=Ws/Wp;
[N,wc]=cheb2ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb2ap(N,rs);
[B,A]=zp2tf(z,p,k);
w=0:0.02*pi:pi;
[h,w]=freqs(B,A,w);
plot(w*wc/wp,20*log10(abs(h)),'k');grid
xlabel('\lambda');
ylabel('A(\lambda)/dB');
%% 椭圆滤波器
% ellipap函数
% [z,p,k]=ellipap(n,Rp,Rs)
% elipord函数
% [n.Wn]=elipord(Wp,Ws,Rp,Rs)
% 返回数字椭圆滤波器的最低阶数和通带截止频率

% ellip函数
% [b,a]=ellip(n,Rp,Rs,Wp)
% 返回具有归一化通带边缘频率Wp的n阶低通数字椭圆滤波器的传递函数系数

%% 7-29 ellipord函数设计椭圆滤波器示例
clear,clc,clf
Wp=3*pi*4*12^3;
Ws=2*pi*12*12^3;
rp=2;
rs=25;          % 设计滤波器的参数
wp=1;ws=Ws/Wp;  % 对参数归一化
[N,wc]=ellipord(wp,ws,rp,rs,'s');   % 计算滤波器阶数和阻带起始频率
[z,p,k]=ellipap(N,rp,rs);  % 计算零点、极点、增益
[B,A]=zp2tf(z,p,k);     %计算系统函数的多项式
w=0:0.03*pi:2*pi;[h,w]=freqs(B,A,w);
plot(w,20*log10(abs(h)),'k');grid 
xlabel('\lambda');ylabel('A(\lambda)/dB');

%% 7-30 ellpap函数设计椭圆滤波器示例
clear,clc,clf
n=0:0.02:4;
for i=1:4
    switch i
        case 1,N=1;
        case 2,N=3;
        case 3,N=5;
        case 4,N=7;
    end
    Rp=1;
    Rs=15;                      %通带波纹1dB，阻带衰减15dB
    [z,p,k]=ellipap(N,Rp,Rs);   %设计椭圆滤波器
    [b,a]=zp2tf(z,p,k);  % 将零/极点增益形式转换为传递函数形式
    [H,w]=freqs(b,a,n)      %按n指定的频率点给出频率响应
    magH2=(abs(H).^2);      %给出传递函数幅度平方
    subplot(2,2,i);plot(w,magH2);grid on
    title(['N=' num2str(N)]);   %将数字N转换为字符串N=
    xlabel('w/wc');
    ylabel(' |H(jw)|^2')
end

% 函数lp2lp用于将模拟低通滤波器转换为实际模拟低通滤波器
% [bt,at]=lp2lp(b,a,wordcloud)
% [At,Bt,Ct,Dt]=lp2lp(A,B,C,D,wordcloud)

%% 7-31 设计合适的切比雪夫I型滤波器，实现低通到低通的频率变换
clear,clc,clf
Wp=3*pi*5000;
Ws=3*pi*13000;
rp=2;
rs=25;
wp=Wp/Wp;ws=Ws/Wp;
[N,wc]=cheb1ord(wp,ws,rp,rs,'s');
[z,p,k]=cheb1ap(N,wc);  
[bp,ap]=zp2tf(z,p,k);
[b,a]=lp2lp(bp,ap,Wp);
w=0:3*pi*120:3*pi*30000;
[h,w]=freqs(b,a,w);     %计算频率响应
plot(w/(2*pi),20*log10(abs(h)),'k');grid
xlabel('f/Hz');
ylabel('幅度/dB');

%% 7-32
clear,clc,clf
Rp=3;
Rs=30;
[z,p,k]=ellipap(4,Rp,Rs);  % 计算零点、极点、增益
[b,a]=zp2tf(z,p,k);
n=0:0.02:4;
[h,w]=freqs(b,a,n);
subplot(121);plot(w,abs(h).^2);grid on
xlabel('w/wc');ylabel(' |H(jw)|^2');
title('原型椭圆滤波器(wc=1)')
[bt,at]=lp2lp(b,a,0.6);
[ht,wt]=freqs(bt,at,n);
subplot(122);plot(wt,abs(ht).^2);grid on
xlabel('w/wc');ylabel(' |H(jw)|^2');
title('原型椭圆滤波器(wc=0.6)')


